The following reduction code works directly from the full set of images in a given nights observations. As observations are added during the evening, the ImageFileCollection object defined at the beginning needs to be refreshed.
Files can be rsync'd from the whtobs computer using the following command:
$ rsync -av whtobs@taurus.ing.iac.es:/obsdata/whta/20170624/ 20170624/
where 20170624/ is set to the relevant observation night.
The starting point for this code was the very helpful scripts of Steve Crawford. However, various changes have been made and key steps (e.g. wavelength calibration/flat processing) differ.
Propagate uncertainties to create variance map
Flux calibration
In [199]:
import numpy as np
from multiprocessing import Pool
from functools import partial
import astropy.units as u
from astropy.io import fits
from astropy.modeling import models, fitting
from astropy.stats import sigma_clip, mad_std
from scipy.ndimage import binary_dilation
from astropy.utils.console import ProgressBar
import ccdproc
from ccdproc import ImageFileCollection, CCDData
def fit_chebyshev(row, degree=5, grow=3):
"""
Fit Chebyshev1D model to a CCD row, including masking of outlier pixels
Params
------
row : array,
Input CCD row to fit polynomial to.
degree : int,
Chebyshev Polynomial Order
grow : int,
Number of iterations to dilate around masked pixels
"""
fitter = fitting.LinearLSQFitter()
input_mask = row.mask
clipped = sigma_clip(row, stdfunc=mad_std)
clipped_pixels = np.array(clipped.mask+row.mask).astype('float')
clipped_pixels = binary_dilation(clipped_pixels, iterations=grow)
row[clipped_pixels==1] = np.median(row)
masked_row = np.ma.array(data=row,
mask=(clipped_pixels == 1),
fill_value=np.median(row))
x = np.arange(len(row))
model = models.Chebyshev1D(degree=degree)
fitted_model = fitter(model, x, row)
return fitted_model(np.arange(len(row)))
def fit_background_old(data, degree=5, grow=3, verbose=True):
"""
Background estimation for longslit CCD image
Params
------
ccd : array,
Input CCD data for background estimation.
degree : int,
Chebyshev Polynomial Order
grow : int,
Number of iterations to dilate around masked pixels
verbose : bool, default=True
If true, print progress bar
"""
fitted_sky = np.zeros_like(data).astype('float')
if verbose:
bar = ProgressBar(data.shape[0], ipython_widget=True)
for irx, row in enumerate(data):
fitted_sky[irx] = fit_chebyshev(row, degree, grow)
if verbose:
bar.update()
return fitted_sky
def fit_background(data, degree=5, grow=3, verbose=True, njobs=4):
"""
Parallelised background estimation for longslit CCD image
Params
------
data : array,
Input CCD data for background estimation.
degree : int,
Chebyshev Polynomial Order
grow : int,
Number of iterations to dilate around masked pixels
njobs : int
Number of processes to initiate for fitting
"""
kwargs={'degree': degree, 'grow': grow}
p = Pool(njobs)
fitted_sky = p.map(partial(fit_chebyshev, **kwargs), data)
p.close()
return np.array(fitted_sky).astype('float')
In [200]:
ic1 = ImageFileCollection('../20170626/')
In [201]:
blue_bias_list = []
for filename in ic1.files_filtered(obstype='Bias', isiarm='Blue arm'):
print ic1.location + filename
ccd = CCDData.read(ic1.location + filename, unit = u.adu)
#ccd = ccdproc.create_deviation(ccd, gain=ccd.header['GAIN']*u.electron/u.adu,
# readnoise=ccd.header['READNOIS']*u.electron)
#this has to be fixed as the bias section does not include the whole section that will be trimmed
ccd = ccdproc.subtract_overscan(ccd, median=True, overscan_axis=0, fits_section='[1:966,4105:4190]')
ccd = ccdproc.trim_image(ccd, fits_section=ccd.header['TRIMSEC'])
blue_bias_list.append(ccd)
master_bias_blue = ccdproc.combine(blue_bias_list, method='median')
master_bias_blue.write('../20170626/processed/master_bias_blue.fits', clobber=True)
red_bias_list = []
for filename in ic1.files_filtered(obstype='Bias', isiarm='Red arm'):
print ic1.location + filename
ccd = CCDData.read(ic1.location + filename, unit = u.adu)
#ccd = ccdproc.create_deviation(ccd, gain=ccd.header['GAIN']*u.electron/u.adu,
# readnoise=ccd.header['READNOIS']*u.electron)
#this has to be fixed as the bias section does not include the whole section that will be trimmed
ccd = ccdproc.subtract_overscan(ccd, median=True, overscan_axis=0, fits_section='[1:966,4105:4190]')
ccd = ccdproc.trim_image(ccd, fits_section=ccd.header['TRIMSEC'] )
red_bias_list.append(ccd)
master_bias_red = ccdproc.combine(red_bias_list, method='median')
master_bias_red.write('../20170626/processed/master_bias_red.fits', clobber=True)
In [203]:
from astropy.convolution import convolve, Gaussian2DKernel
kernel = Gaussian2DKernel(25)
red_flat_list = []
for filename in ic1.files_filtered(obstype='Flat', isiarm='Red arm', object='good flat red'):
ccd = CCDData.read(ic1.location + filename, unit = u.adu)
#ccd = ccdproc.create_deviation(ccd, gain=ccd.header['GAIN']*u.electron/u.adu,
# readnoise=ccd.header['READNOIS']*u.electron)
#this has to be fixed as the bias section does not include the whole section that will be trimmed
ccd = ccdproc.subtract_overscan(ccd, median=True, overscan_axis=0, fits_section='[1:966,4105:4190]')
ccd = ccdproc.trim_image(ccd, fits_section=ccd.header['TRIMSEC'] )
ccd = ccdproc.subtract_bias(ccd, master_bias_red)
red_flat_list.append(ccd)
master_flat_red = ccdproc.combine(red_flat_list, method='median')
convolved_flat_red = convolve(master_flat_red.data, kernel, boundary='extend')
master_flat_red.write('../20170626/processed/master_flat_red.fits', clobber=True)
master_flat_red.data /= convolved_flat_red
master_flat_red.write('../20170626/processed/master_flat_norm_red.fits', clobber=True)
In [204]:
blue_flat_list = []
for filename in ic1.files_filtered(obstype='Flat', isiarm='Blue arm', object='good flat blue'):
ccd = CCDData.read(ic1.location + filename, unit = u.adu)
#ccd = ccdproc.create_deviation(ccd, gain=ccd.header['GAIN']*u.electron/u.adu,
# readnoise=ccd.header['READNOIS']*u.electron)
#this has to be fixed as the bias section does not include the whole section that will be trimmed
ccd = ccdproc.subtract_overscan(ccd, median=True, overscan_axis=0, fits_section='[1:966,4105:4190]')
ccd = ccdproc.trim_image(ccd, fits_section=ccd.header['TRIMSEC'] )
ccd = ccdproc.subtract_bias(ccd, master_bias_blue)
blue_flat_list.append(ccd)
master_flat_blue = ccdproc.combine(blue_flat_list, method='median')
convolved_flat_blue = convolve(master_flat_blue.data, kernel, boundary='extend')
master_flat_blue.write('../20170626/master_flat_blue.fits', clobber=True)
master_flat_blue.data /= convolved_flat_blue
master_flat_blue.write('../20170626/processed/master_flat_norm_blue.fits', clobber=True)
In [314]:
ic1 = ImageFileCollection('../20170626/')
objects = ['GRG3'] #'SP1234+254', 'USS439'] 'USS268'
for objname in objects:
print(objname)
"""
Blue Arm
"""
blue_target_list = []
for ifx, filename in enumerate(ic1.files_filtered(obstype='TARGET', isiarm='Blue arm', object=objname)):
print('{0} {1}'.format(ifx+1, filename))
hdu = fits.open(ic1.location + filename)
ccd = CCDData(hdu[1].data, header=hdu[0].header+hdu[1].header, unit = u.adu)
#this has to be fixed as the bias section does not include the whole section that will be trimmed
ccd = ccdproc.subtract_overscan(ccd, median=True, overscan_axis=0, fits_section='[1:966,4105:4190]')
ccd = ccdproc.trim_image(ccd, fits_section=ccd.header['TRIMSEC'] )
ccd = ccdproc.subtract_bias(ccd, master_bias_blue)
ccd = ccdproc.flat_correct(ccd, master_flat_blue)
ccd = ccdproc.cosmicray_lacosmic(ccd, sigclip=4.5, niter=10, sigfrac=0.05, psffwhm=2.5,
gain=ccd.header['GAIN'], readnoise=ccd.header['READNOIS'])
# Do sky subtraction
ccd.mask[:,785:800] = True
sky = fit_background(np.ma.array(ccd.data, mask=ccd.mask))
ccd.data -= sky
# Rotate Frame
ccd.data = ccd.data.T
ccd.mask = ccd.mask.T
blue_target_list.append(ccd)
#ccd.write('obj_'+filename, clobber=True)
blue_target = ccdproc.combine(blue_target_list, method='average')
blue_target.write('../20170626/processed/{0}_blue.fits'.format(blue_target_list[0].header['object']), clobber=True)
for objname in objects:
print(objname)
"""
Red Arm
"""
red_target_list = []
for ifx, filename in enumerate(ic1.files_filtered(obstype='TARGET', isiarm='Red arm', object=objname)):
print('{0} {1}'.format(ifx+1, filename))
hdu = fits.open(ic1.location + filename)
ccd = CCDData(hdu[1].data, header=hdu[0].header+hdu[1].header, unit = u.adu)
#this has to be fixed as the bias section does not include the whole section that will be trimmed
ccd = ccdproc.subtract_overscan(ccd, median=True, overscan_axis=0, fits_section='[1:966,4105:4190]')
ccd = ccdproc.trim_image(ccd, fits_section=ccd.header['TRIMSEC'] )
ccd = ccdproc.subtract_bias(ccd, master_bias_red)
ccd = ccdproc.flat_correct(ccd, master_flat_red)
ccd = ccdproc.cosmicray_lacosmic(ccd, sigclip=4.5, niter=10, sigfrac=0.05, psffwhm=2.5,
gain=ccd.header['GAIN'], readnoise=ccd.header['READNOIS'])
# Do sky subtraction
ccd.mask[:,690:700] = True
sky = fit_background(np.ma.array(ccd.data, mask=ccd.mask))
ccd.data -= sky
# Rotate Frame
ccd.data = ccd.data.T
ccd.mask = ccd.mask.T
#ccd.write('obj_'+filename, clobber=True)
red_target_list.append(ccd)
red_target = ccdproc.combine(red_target_list, method='average')
red_target.write('{0}_red.fits'.format(red_target_list[0].header['object']), clobber=True)
red_target.mask[785:800,:] = True
red_sky = fit_background(np.ma.array(red_target.data.T, mask=red_target.mask.T)).T
red_target.data -= red_sky
red_target.write('../20170626/processed/{0}_red.fits'.format(red_target_list[0].header['object']), clobber=True)
In [38]:
for filename in ic1.files_filtered(obstype='Arc', isiarm='Blue arm', object='CuNe+CuAr b tar'):
hdu = fits.open(ic1.location + filename)
ccd = CCDData(hdu[1].data, header=hdu[0].header+hdu[1].header, unit = u.adu)
#this has to be fixed as the bias section does not include the whole section that will be trimmed
ccd = ccdproc.subtract_overscan(ccd, median=True, overscan_axis=0, fits_section='[1:966,4105:4190]')
ccd = ccdproc.trim_image(ccd, fits_section=ccd.header['TRIMSEC'] )
ccd = ccdproc.subtract_bias(ccd, master_bias_blue)
ccd = ccdproc.flat_correct(ccd, master_flat_blue)
ccd = ccdproc.cosmicray_lacosmic(ccd, sigclip=4.5, gain=ccd.header['GAIN'], readnoise=ccd.header['READNOIS'])
ccd.data = ccd.data.T
ccd.mask = ccd.mask.T
ccd.write('arc_blue_'+filename, clobber=True)
for filename in ic1.files_filtered(obstype='Arc', isiarm='Red arm', object='CuNe+CuAr r tar'):
hdu = fits.open(ic1.location + filename)
ccd = CCDData(hdu[1].data, header=hdu[0].header+hdu[1].header, unit = u.adu)
#this has to be fixed as the bias section does not include the whole section that will be trimmed
ccd = ccdproc.subtract_overscan(ccd, median=True, overscan_axis=0, fits_section='[1:966,4105:4190]')
ccd = ccdproc.trim_image(ccd, fits_section=ccd.header['TRIMSEC'] )
ccd = ccdproc.subtract_bias(ccd, master_bias_red)
ccd = ccdproc.flat_correct(ccd, master_flat_red)
ccd = ccdproc.cosmicray_lacosmic(ccd, sigclip=4.5, gain=ccd.header['GAIN'], readnoise=ccd.header['READNOIS'])
ccd.data = ccd.data.T
ccd.mask = ccd.mask.T
ccd.write('arc_red_'+filename, clobber=True)
This section probably needs to be set up so that the wavelength calibration is done for the image section extracted for the target spectrum - with the calibration done from the appropriate arcs. i.e. load in the correct corresponding arcs based on position/time of observation and calculate the wavelength solution.
For now, this just calculates a rough solution based on one example.
In [127]:
Fig, Ax = plt.subplots(2, 1, figsize=(12,8))
blue_arc_1d = np.median(CCDData.read('arc_blue_r2545050.fit').data/convolved_flat_blue.T, axis=0)
Ax[0].plot(np.arange(4099), blue_arc_1d)
Ax[0].set_xlim([1000,4000])
#Ax[0].set_ylim([0, 2000])
red_arc_1d = np.median(CCDData.read('arc_red_r2545051.fit').data/convolved_flat_red.T, axis=0)
Ax[1].plot(np.arange(4096), red_arc_1d)
Ax[1].set_xlim([1000,4000])
#Ax[1].set_ylim([0, 30000])
Fig.tight_layout()
# Identified lines for fitting and their corresponding wavelengths.
# These may not be accurate enough alone, so the plan is to use these to provide the starting
#
blue_lines_pix = [1177, 1868, 1895, # 1846,
2140, 2173, 2368, 2610,
2703, 2950, 3000, 3132,
3556, 3638]
blue_lines_w = [3273.96, 3868.53, 3891.98, # 3858.58,
4103.91, 4131.72, 4300.1, 4510.73,
4589.9, 4806.02, 4847.81, 4965.08,
5330.78, 5400.56]
red_lines_pix = [1341, 1392, 1938, 1994,
2057, 2072, 2111, 2188,
2498, 2672, 2814, 3143,
3199, 3270, 3437, 3505]
red_lines_w = [5852.49, 5944.83, 6929.47, 7032.41,
7147.04, 7173.94, 7245.17, 7383.98,
7948.17, 8264.52, 8521.44, 9122.97,
9224.5, 9354.22, 9657.78, 9784.5]
In [177]:
def line_models(mean_init):
lmodels = []
for ilx, line in enumerate(mean_init):
compound_model = models.Gaussian1D(mean=mean_init[ilx],
amplitude=1.,
stddev=1.,
bounds = {'stddev':(0,5)})
lmodels.append(compound_model)
return lmodels
def fit_lines(spectra, mean_init):
""" Fit emission lines for an observed WL target
"""
lmodels = line_models(mean_init)
fit = fitting.LevMarLSQFitter()
fitted_models = []
amplitudes = []
means = []
widths = []
X = np.arange(len(spectra))
#plt.plot(ccd_wavelengths.value, obs_flux[i])
for model in lmodels:
#dprint(wcenter)
try:
lower = int(np.maximum(0, model.mean-10))
upper = int(np.minimum(len(spectra), model.mean+10))
fitted_model = fit(model, X[lower:upper], spectra[lower:upper])
fitted_models.append(fitted_model)
amplitudes.append(fitted_model.amplitude.value)
means.append(fitted_model.mean.value)
widths.append(fitted_model.stddev.value)
#plt.plot(ccd_wavelengths.value[lower:upper], fitted_model(ccd_wavelengths.value[lower:upper]))
except:
fitted_models.append(model)
amplitudes.append(-99.)
means.append(model.mean.value)
widths.append(model.stddev.value)
raise
return means, widths, amplitudes, fitted_models
def evaluate_models(X, models):
model_flux = np.zeros_like(X)
for model in models:
model_flux += model(X)
return model_flux
def get_solution(spectrum, reference_lines_pix, reference_lines_lambda):
means, widths, amplitudes, fmodels = fit_lines(spectrum, reference_lines_pix)
fit = fitting.LinearLSQFitter()
m = models.Legendre1D(degree=3)
solution = fit(m, means, reference_lines_lambda)
solution_lambda = solution(means)
pixels = np.arange(len(spectrum))
rms = np.sqrt(np.mean((solution_lambda-reference_lines_lambda)**2))
return solution(pixels), rms, solution
def wavelength_calibration(image, reference_lines_pix, reference_lines_lambda,
ref_file1, ref_file2):
"""
Initial Solution
"""
wavelengths = np.zeros(image.shape)
rms = np.zeros(image.shape[0])
for irx, row in enumerate(image):
wavelengths[irx], rms[irx], = get_solution(row, reference_lines_pix, reference_lines_lambda)
"""
Second solution
in progress
"""
return wavelengths, rms
In [179]:
data = CCDData.read('arc_blue_r2545050.fit')
data.data = test
data.write('test_wavelength_solution.fits', overwrite=True)
In [136]:
means, widths, amplitudes, fmodels = fit_lines(blue_arc_1d, blue_lines_pix)
X = np.arange(len(blue_arc_1d)).astype('float')
mf = evaluate_models(X, fmodels)
fit = fitting.LinearLSQFitter()
m = models.Legendre1D(degree=3)
fwave = fit(m, means, blue_lines_w)
%matplotlib notebook
plt.plot(means, blue_lines_w, 'o')
plt.plot(X, fwave(X))
Out[136]:
In [137]:
means, widths, amplitudes, fmodels = fit_lines(red_arc_1d, red_lines_pix)
X = np.arange(len(red_arc_1d)).astype('float')
mf = evaluate_models(X, fmodels)
fit = fitting.LinearLSQFitter()
m = models.Legendre1D(degree=3)
fwave = fit(m, means, red_lines_w)
%matplotlib notebook
plt.plot(means, red_lines_w, 'o')
plt.plot(X, fwave(X))
Out[137]:
In [313]:
from astropy.convolution import Gaussian1DKernel
# Wavelength Solution
data_blue = CCDData.read('arc_blue_r2545050.fit')
test_blue, rms = get_solution(data_blue.data.sum(0), blue_lines_pix, blue_lines_w)
data_red = CCDData.read('arc_red_r2545051.fit')
test_red, rms = get_solution(data_red.data.sum(0), red_lines_pix, red_lines_w)
# Smoothing kernel, if required.
kernel = Gaussian1DKernel(stddev=3)
data_blue = CCDData.read('../20170626/processed/GRG3_blue.fits')
# Range needs changed to match the region of interest as it can shift slightly for different objects
# GRGs are more extended so it can be larger. The offset applied to the _sky array is arbitrary currently.
blue_1d = data_blue.data[504:534,:].sum(axis=0)
blue_smoothed = convolve(blue_1d, kernel)
blue_smoothed_sky = convolve(data_blue.data[502+30:529+30,:].sum(axis=0), kernel)
data_red = CCDData.read('../20170626/processed/GRG3_red.fits')
red_1d = data_red.data[412:465,:].sum(axis=0)
red_smoothed = convolve(red_1d, kernel)
red_smoothed_sky = convolve(data_red.data[412-100:465-100,:].sum(axis=0), kernel)
Fig, Ax = plt.subplots(1)
#Ax.plot(np.arange(len(blue_smoothed)), blue_smoothed-blue_smoothed_sky)
#Ax.plot(np.arange(len(red_smoothed)), red_smoothed-red_smoothed_sky)
Ax.plot(test_blue, blue_smoothed) #-blue_smoothed_sky
Ax.plot(test_red, red_smoothed) #-red_smoothed_sky)
Out[313]: